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6-SEP=1984 11 MTYRTL. SRCIMTHDSORT.MAR; 1 (1) 
-TITLE MTHSDSQRT 3 Double, Floating Point Square Root routine 
. IDENT /1-015/ : File: MTHDSQRT.MAR EDIT: RNH1015 
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COPYRIGHT (c) 1978, 1980, 1982, 1984 BY 
DIGITAL EQUIPMENT CORPORATION, MAYNARD, MASSACHUSETTS. 
ALL RIGHTS RESERVED. 


THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED 
: ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE 
Pp 


s 
ca cy 
% a 
* » 
oa oa 
e 2 
* 8 
* ON OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER * 
* THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY” * 
* ERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY * 
* TRANSFERRED. : 
& a 
* THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOT 3 * 
* ® 
® 
e ey 
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* e 
& ry 
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IC 
AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMEN 
CORPORATION. 


DIGITAL ASSUMES NO RESPONSIBI 
SOFTWARE ON EQUIPMENT WHICH I 


LITY FOR THE USE OR RELIABILITY OF ITS 
S NOT SUPPLIED BY DIGITAL. 
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FACILITY: MATH LIBRARY 
++ 
ABSTRACT: 


MTHSDSQRT is a function which returns the floating point square root 
of its single precision floating point argument. The call is standard 
TAN i] oy rence. : ‘ 

MTHSDSQRT_RS is a special routine which is the same as MTH$DSQRT 
except a faster non-standard JSB call is used with the argument in 
RO and no registers are saved. 


VERSION: 01 
HISTORY: 
AUTHOR: 
Peter Yuo, 15-Oct-76: Version 01 
MODIFIED BY: 


Ql-1 Peter Yuo, gf -nay- 77 
01-2 Peter Yuo, 31-May-77 
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et ted barrent Edit History at et 91:35:98 EMTHRIL. SREIMTHDSORT MAR: 1 . (3) 


~SBTTL HISTORY ; Detailed Current Edit History 


MTHSDSQRT 
1-015 


“eo 
=o 


oc 
aoe 


x-- 
i om 
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ALGORITHMIC DIFFERENCES FROM FP=11/C ROUTINE: 
1. Last DIVD is rounees instead of truncated. Results should be 
correct within 2 LSB's. 


Edit History for Version 01 of MTHSDSQRT 


gin) Code saving after code review March 1977 
1-2 ROTL shift in gerbone into highest bit. Use ASHL instead. 
ADDL instruction after ADJUST has begn changed into ADDW to prevent 
overflow if R1<31:16> = FFFF and RO<31:16> = FFFF 


01-3 Finish error handling 10-June-1977 
0-4 MTHSSERROR anaes to MTH$$SIGNAL. 
MTH$_... changed to MTH__.... 
Changed error handling mechanism. Put error result in RO:R1 before 
calling MTHSS$SIGNAL in order to allow user modify error result. 
Return -0.0 on negative arg. TNH 20-Dec-77 
Remove undefined $1 lobal. TNH i pt, 
d Rich Lary's code bums. JMT 26-Jan-78 
Move .ENTRY symbol definition to module header. TNH 14-Aug-78 
peeete version number and copyrtane notice. JBS 16-NOV-78 
Change symbol MTH_ SQUROONEG to PMTHSK SQUROONEG. JBS 07-DEC-78 
Add "_"' to the PSECT directive. JBS 22-DEC-78 
Declare externals. SBL 17-May-1979 
Use ASHQ, mot ASHL, to generate reserved operand. JAW 12-Oct-1979 
Changed W* to G* in call to MTH$$SIGNAL RNH 09-Sept-1981 
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3 beciaraties Part of Modul ae oets 91:33:98 f at e° 


-SBTTL DECLARATIONS ; Declarative Part of Module 


; INCLUDE FILES: MTHJACKET.MAR 


: EXTERNAL SYMBOLS: 


-DSABL GBL 
-EXTRN MTHSK_SQUROONEG 
eEXTRN MTHSSSIGNAL ; SIGNAL SEVERE error 


; EQUATED SYMBOLS: 


ACMASK = “M<IV, R2, R3, R4, R5> ; register save mask and IV enable 


MACROS: none 
PSECT DECLARATIONS: 
-PSECT _MTHSCODE PIC,SHR,LONG,EXE,NOWRT 


3 program section for math routines 


OWN STORAGE: none 
CONSTANTS: 


: Constants A and B chosen for k = odd 


LF_ODD_A_E63 “X13CD5FD4 
LFODD~B7EM63 “%3C4A2018 


Constants A and B chosen for k = even 


“XF61A4015 
“X4B231FD7 


LF_EVEN_A 
LFEVENTB_EM64 


MTHRTL.SRCIMTHDSORT.MAR; 1 
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Mors ondard Double Precision F gree b= 984 9:33:98 YATHRTL. SRE TATNDSORT «MAR: 1 ‘soe (a) 


~-SBTTL MTHSDSQRT = Standard Double Precision Floating DSQRT 


le 
QRT 


{FUNCTIONAL DESCRIPTION: 

; DSQRT = double precision floating point function 

DSQRT(X) is computed using the following approximation technique: 
If X <= 0, error. Let X = IXi. 
Let X = 2**K * F where F is the fractional part. 


If K = even, X = 2e2(2P) * F, 
DSQRT(X) = 2e*P * DSORT(F), 1/2 =< F <1 


If K = odd, X = 2e#(2P+1) # F = 2ee(2P+2) * (F/2), 
DSQRT(X) = 2ee(P4+1) ® DSQRT(F/2), 1/4 =< F/2 < 1/2. 


Let F° = AtF + B, 


WANIANAIAINIIIrononn) Ue 


OOoOooooooooooooooooooo 


octal), 
octal), for K = even. 


or 


At(F/2) + B 


octal), 
octal), for K = odd. 


ow el 


or 


and 


K' for K = 


oP, 
=P +1 for K 
Let YO = 2**k' * F* as a ight Line eogres iaetton wthin the 
given interval using coefficients A and B which minimize the 
absolute error at the midpoint and endpoint. 
Starting with YO, three Newton-Raphson iterations are performed. 
YCn+1) = (1/2) * ( YEn) + xX/¥Efn]) 
The relative error is < 10*#*-17, 

CALLING SEQUENCE: 

dsaqrt.wd.v = MTHSDSQRT(x.rd.r) 

INPUT PARAMETERS: 


00000004 


00000004 ; define longword multiplier 


LONG = 4 w 

x = 1 * LONG ; Contents of x is the argument 
IMPLICIT INPUTS: none 
OUTPUT PARAMETERS: 

VALUE: double precision square root of the argument 
IMPLICIT OUTPUTS: none 


MRO OONAUES WIN (OOD NAUES WN OOO NAU EWN O OO NAME WN O OONOUS WN OOOnNIG 
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MTHSDSOQRT ; Double Floating Point Square Root rout 16-SEP-1984 01:22: AX/VMS Macro v04-00 Pa 5 
vat ATHSDSORT - Standard Double Precision F mie} 7 97:53:96 EMTHRTL. SREJMTHDSORY MAR: 1 - (4) 
00 ! ; ; COMPLETION CODES: none 
09 5 ; SIDE EFFECTS: 
00 1 5 : Signals: MTH$_SQUROONEG if X < 0.0 with reserved operand in RO/R1 
4 1 8 ; (copied to the signal mechanism vector CHF$L_MCH_RO/R1 by LIBSSIGNAL). 
0 189 ; Associated message is: ‘SQUARE ROOT OF NEGATIVE VALUE". Result is reserved 
8 139 ; operand -0.0 unless a user supplied (or any) error handler changes CHF$L_MCH_RO/R1 
0000 136 : NOTE: This procedure disables floating point underflow, enables integer 
000 1935 ; overflow, causes no floating overflow or other arithmetic traps, and 
$888 194 ; preserves enables across the call. 
0 195 ; 
0000 196 ;--- 
0000 197 
0000 198 
403C 0000 199 «ENTRY MTHSDSQRT, ACMASK 3; standard call-by-reference entry 
BooS 200 ; disable DV (and FU), enable Iv 
B28 201 MTHSFLAG_JACKET ; flag that this is a jacket procedure in 
6D 00000000 ' GF 9E $095 MOVAB G*MTHSSJACKET_HND, (FP) 
0009 3; set handler address to jacket 
8 3; handler 
0 f ; case of an error in special routine 
0 0 MOVD ax (AP), RO ; RO/R1 = org 
0 04 BSBB MTHSDSQRT_RS ; call kernel DSQRT rountine 
0 05 RET 3 return - result in RO/R1 
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| MTHSDSQRT ; Double Floating Point Square Root rout 16-SEP-1984 01:22: AX/VMS Macro V04-00 Pa eT 
8 ATHSDSORT_RS - Special SORT routine oreeen19bs 91:55:50 EMTURTE. SaeamrHDsanT.mar:1 § 2% (8 Le 
! 4 -SBTTL MTHSDSQRT_RS = Special DSQRT routine 
: 09 : Special DSQRT - used by the standard routine, and directly. 
01 11 $ CALLING SEQUENCE: 
01 \¢ ; save anything in R2:R5 . 
4 1s 3 0 saan © : input in RO/R1 
1 14; JSB MTHSOSQRT_RS 
001 15 ; return with resul€ in RO/R1 
001 16 ; Note: This routine is written to avoid causing any integer overflows, floating 
091 18 3 suet varts floating underflows or divide by 0 conditions, whether enabled or 
3 not. 
1 19 ; 
0010 0 ; REGISTERS USED: 
001 1 3 RO/R1 = Floating argument then result 
001 § : R2/R3 = scratc 
Bob 7 3 R4/R5 = hold X during calc of F*, K*. 
$010 § 5 MTHSDSQRT_RS:: ; JSB routine for DSQRT 
54 50 70 0010 6 mOvD RO, R4 3 test sign of X and save it in R4/R5. 
5p Ss «15 Bote st h BLEQ ZERO_NEG ; branch to ZERO_NEG if X =< 0 
ge gee 
51 50 3¢ 0015 31 Pos: MOVZWL RO, R1 ; isolate low 16 bits (sign,exp,>fract) in R 
51 94 0018 $35 CLRB R1 3; R1 now has sign and left 7 exp bits 
50 51 AA OO1A 3 BICW Ri, RO ; clear sign and left 7 exp bits 
50 95 Q01D 234 TSTB RO : check low bit of exp 
10 18 OO1F $32 BGEQ EVEN : and branch if 1 
50 i3CD5FD4 BF 44 0021 36 MULF #LF_ODD_A_E63, RO ; add 64 (half of bias) to (exponent-2) 
0028 237 3 and start approximation calc 
50 3C€4A2018 8F 40 0028 238 ADDF #LF_ODD_B_EM63, RO ; RO = (first approx) * 2**-64 
aw poet $73 ERB ADJOST 3 go adjust 
0031 40 
bos) 41 EVEN: ; 
50 2000 8F AO 0031 $s ADDW #*¥°2000, RO : exp is 0 - make it 64 (2**-64) for legalit 
50 F61A4015 BF 44 0036 4 MULF #L: =VEN_A, RO 
50 4B231FD7 BF 40 Fi rf: snsuee ADDF #LF_EVEN_B_EM64, RO ; RO = (first approx) * 2**-64 
51 51 1F 9C 0044 $2 " ROTL #31, R1, R1 ; divide R1 (exptbias) by 2, 
st 4 ste : giving (exp/2+64) 
50 51 AO 00468 248 ADDW = Ri, RO ; Insert exp/2 in first approx and 
Bhee rs 3; re-bias it. 
i 2) ; first iteration - single precision is sufficient 
51 54 50 47 0048 34 : DIVF3 RO, R4, RI : Rl = 
50 51 40 O04F 54 ADDF . 3; RO = YO + x/Y0 
50 0080 8F A2 O36 55 SUBW #*x80, RO : RO = Y1 = (1/2)(Y0 + x/Y¥0) 
p25 28 : no overflow possible 
057 28 3 second iteration, do in double precision to get truncated( rather than 
$2) 23 : rounded) result. 
51 D4 08 61° CLRL_ RI : lower part (x) = 0 
52 54 50 67 0059 6¢ DIVD3 RQ. R4, R2 : ace = x/Y1 
50 52 60 005D 6 ADDD R2, RO ; RO/R1 = Y1 + higher_part(Xx/Y1) 


i 


7E OO'8F 9A 


50 01 OF 79 condition value 


RO/R1 = result = reserved operand -0.0 
RO/R1 goes to si nal mechanism vector 
CCHFSL MCH_RO/R1) so error handler 

con modify the BF, 


PUSHL (SP) 
MOVZBL #MTHS$K_SQUROONEG, -(SP) 
ASHQ #15, #T, RO 


OCOO0000'GF 02 FB CALLS #2, G*MTHSSSIGNAL 
RSB 
END 


oe error and use real user's PC 
independent of CALL vs JSB 
return - RO pestered. tren CHF SL_MCH_RO/R1 


Bete Ge Ge Ge Se Se Ge Se 


MTHSDSOQRT ; Double Floating Point Square ies font bp AX/VMS Macro V04-00 Pa 7 
arate ATHSDSORT_RS. = Speci ial DSQRT Bran + Bs at 9 9 ‘$3: 96 EMTHRIL. SREIMTMDSORT MAR: 1 ” (5 
50 0080 8F A2 p 3 rt: SUBW #*x80, RO 3; RO/R1 = Y2 = (1/2) (Y14#xX/Y¥1) 
6 ; third iteration, do in double precision. 
> a | 29 67 be3 68 Divbd3 RQ. R4, R2 3: R2/R3 = X/Y 
50. 5 6 069 $? ADDD R2, RO ; 4 = v2+ X/Y 
50 0080 8F A Bes 0 SUBW #*x80, RO 3; RO/R1 = DSQRT(X) = 
8 1 71 3 (172) (¥24x/Y2) 
05 0071 i SQRTX: RSB 3 return, RO/R1 = result 
007 7 
007 74; % =< 0 
07 75 ; 
07 i: ZERO_NEG: 
FD 13 007 7 BEQL SQRTX ; return with RO = result = 0 
6—E DD 0074 78 return PC from JSB routine 
0076 79 
007A 80 
007E 81 
007E 6 
007E 8 
443 84 
008 85 
0085 86 
0086 287 
86 288 
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| MTHSDSQRT ; Double Floating Point Square Root rout bp a BD 01:58:06 ae Macro V04-00 Page § 
Symbol table 6-SEP-1984 11:22: MTHRTL.SRCJMTHDSORT MAR; 1 (5) 
ACMASK = 4 ed tg 
ADJUST 00000044 R 4 
N sere eth R 1 
LF _EVEN_A = FOIAG 13 
LF_EVEN_B EM64 = sReaire 
LF-ODD_A_E63_ = 13CDSFD4 
LF ODD_B_EM63 = 3C€4A2018 
LONG = 00000004 
MIHSSJACKET_HND teeeeeee x 01 
MTHSSSIGNAL ereerere X 00 
MTHSDSQRT 00000000 RG 01 
MTHSDSORT RS 00000010 RG Ba 
MTHSK_SQUROONEG teeneeee x 0 
POS 00000015 R 01 
SQRTX 00000071 R 01 
x = 00000004 
ZERO_NEG 00000072 R 01 
os + 
: Psect synopsis ! 
to nwe rma woe ren scoe + 
PSECT name Allocation PSECT No. Attributes 
. ABS, 00000000 «¢ 0.) 00 ¢ 0.) NOPIC USR CON ABS LCL NOSHR NOEXE NORD NOWRT NOVEC BYTE 
_MTHSCODE 00000086 (¢ 134.) O01 ¢ 1.) PIC USR CON REL LCL SHR EXE RD NOWRT NOVEC LONG 
Pees reer ewe nemo naman em eee + 
; Performance indicators ! 
Phase Page faults CPU Time Elapsed Time 
Initialization 9 00: 90:00. 08 00:00:01.27 
Command processing 106 00:00:00.6 00:00:05.92 
Pass 82 00:00:00.88 00:00:04.28 
Symbol table sort eS SE 00:00:00.01 
Pass 2 64 00: 8708-88 00:00:03.11 
Symbol table output 3 00:00: 3-2 00:00:00.16 
Psect synopsis output 3 00:60:00.01 00:00:00.02 
Cross-reference output as So 00:00:00.00 
Assembler run totals 288 00:00:02. 00:00:14.77 
~ working set Limit was 900 —_ 
4182 byies (9 pages) of virtual memory were used to buffer the intermediate code 


& 
There were 10 pages of symbol table space allocated to hold 17 non-local and 0 local symbols. 
348 source Lines were read in Pass 1, producing 11 object records in Pass 2. 

1 page of virtual memory was used to define 1 macro. 


Sere weer eee seer es eer erm en ae} 


Macro Library name Macros defined 


~-$255$DUA28: CSYSLIBISTARLET.MLB; 2 0 


; Double Floating Point Square Rost aut ba 35 ool =138¢ 91: 5: 9¢ ales Ge o V04-00 Page (2) 


aT QRT 
warts Macro Run Statistics MTHRTL SRE SMTHDSORT MAR: 1 


0 GETS were required to define 0 macros. 
There were no errors, warnings or information messages. 
MACRO/ENABLE=SUPPRESS1ON/D1ISABLE=(GLOBAL , TRACEBACK) /LIS=LIS$:M7SDSQRT/OBJ=OBJ$:MTHDSQRT MSRC$:MTHJACKET/UPDATE=(ENHS :MTHJACKET) +MSRC 
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